home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / ARGONET / PD / MATHS / RLAB / RLAB125.ZIP / !RLaB / toolbox / corr < prev    next >
Text File  |  1995-02-19  |  2KB  |  78 lines
  1. //-------------------------------------------------------------------//
  2. // Synopsis:    Compute correlations faster than DSP toolbox, 
  3. //              and with bandpassing
  4.  
  5. // Usage:
  6.  
  7. //   y = corr(x)
  8. //   Compute the autocorrelation of vector x.  Does NOT handle
  9. //   matrices  the way xcorr (in the Signal Processing toolbox) does;
  10. //   instead, a matrix is treated as a vector (with reshape). 
  11.  
  12. //   y = corr(x1,x2)
  13. //   Return the cross-correlation of vectors x1 and x2, which need not
  14. //   be  the same length.  This is similar to xcorr, but faster (the
  15. //   latter efficiently computes 4 correlations and throws 3 away).
  16. //   Also, if  length(x1) is different from length(x2), xcorr returns
  17. //   a wrong-length  result.
  18.  
  19. //   y = corr(x1,x2,f) 
  20. //   As above, but also filter by the bandpass filter f.  f is a
  21. //   two-element vector, with values in [0,1], that specifies the
  22. //   passband; in f, the value 1 represents the Nyquist frequency.
  23.  
  24. // 2/94 David K. Mellinger  dkm1@cornell.edu
  25. // This file is a translation of corr.m from the Osprey toolbox.
  26.  
  27. require isreal
  28. require nextpow2 
  29.  
  30. //-------------------------------------------------------------------//
  31.  
  32. corr = function (in1, in2, f)
  33. {
  34.   local (in1, in2, f)
  35.  
  36.   if (!exist (in2)) { in2 = in1; }
  37.  
  38.   x1 = in1[:];
  39.   x2 = in2[:];
  40.   n1 = length (x1);
  41.   n2 = length (x2);
  42.  
  43.   // pad with 0's for fft & circular corr
  44.   nfft = 2 * 2^nextpow2 (max (n1,n2));
  45.  
  46.   y = conj (fft ([x1;zeros(nfft-n1,1)])) .* fft([x2;zeros(nfft-n2,1)]);
  47.   if (exist (f))
  48.   {
  49.     if (length(f) == 1)
  50.     {
  51.       f = [0, f];    // one number means lowpass filter
  52.     }
  53.     n = nfft/2;
  54.     i0 = round(f[1] * n);
  55.     i1 = round(f[2] * n);
  56.     y[   1 : i0] = zeros (i0,1);
  57.     y[i1+1 : n]  = zeros (n-i1,1);
  58.     y[nfft+1-i0 : nfft]    = zeros (i0,1);
  59.     y[nfft+1-n  : nfft-i1] = zeros (n-i1, 1);
  60.   }
  61.  
  62.   y = ifft (y);
  63.   y = [y[nfft+2-n1:nfft]; y[1:n2]];
  64.   
  65.   // Make output have same form as input.
  66.   if (all (isreal (in1)) && all (isreal (in2)))
  67.   {
  68.     y = real(y);
  69.   }
  70.  
  71.   if (in1.nr == 1)
  72.   {
  73.     y = y.'; 
  74.   }
  75.  
  76.   return y;
  77. };
  78.